home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / surd.cal < prev    next >
Text File  |  1995-07-17  |  4KB  |  289 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Calculate using quadratic surds of the form: a + b * sqrt(D).
  7.  */
  8.  
  9. obj surd {a, b};        /* definition of the surd object */
  10.  
  11. global surd_type = -1;        /* type of surd (value of D) */
  12. static obj surd surd__;        /* example surd for testing against */
  13.  
  14.  
  15. define surd(a,b)
  16. {
  17.     local x;
  18.  
  19.     obj surd x;
  20.     x.a = a;
  21.     x.b = b;
  22.     return x;
  23. }
  24.  
  25.  
  26. define surd_print(a)
  27. {
  28.     print "surd(" : a.a : ", " : a.b : ")" :;
  29. }
  30.  
  31.  
  32. define surd_conj(a)
  33. {
  34.     local    x;
  35.  
  36.     obj surd x;
  37.     x.a = a.a;
  38.     x.b = -a.b;
  39.     return x;
  40. }
  41.  
  42.  
  43. define surd_norm(a)
  44. {
  45.     return a.a^2 + abs(surd_type) * a.b^2;
  46. }
  47.  
  48.  
  49. define surd_value(a, xepsilon)
  50. {
  51.     local    epsilon;
  52.  
  53.     epsilon = xepsilon;
  54.     if (isnull(epsilon))
  55.         epsilon = epsilon();
  56.     return a.a + a.b * sqrt(surd_type, epsilon);
  57. }
  58.  
  59. define surd_add(a, b)
  60. {
  61.     local obj surd    x;
  62.  
  63.     if (!istype(b, x)) {
  64.         x.a = a.a + b;
  65.         x.b = a.b;
  66.         return x;
  67.     }
  68.     if (!istype(a, x)) {
  69.         x.a = a + b.a;
  70.         x.b = b.b;
  71.         return x;
  72.     }
  73.     x.a = a.a + b.a;
  74.     x.b = a.b + b.b;
  75.     if (x.b)
  76.         return x;
  77.     return x.a;
  78. }
  79.  
  80.  
  81. define surd_sub(a, b)
  82. {
  83.     local obj surd    x;
  84.  
  85.     if (!istype(b, x)) {
  86.         x.a = a.a - b;
  87.         x.b = a.b;
  88.         return x;
  89.     }
  90.     if (!istype(a, x)) {
  91.         x.a = a - b.a;
  92.         x.b = -b.b;
  93.         return x;
  94.     }
  95.     x.a = a.a - b.a;
  96.     x.b = a.b - b.b;
  97.     if (x.b)
  98.         return x;
  99.     return x.a;
  100. }
  101.  
  102.  
  103. define surd_inc(a)
  104. {
  105.     local    x;
  106.  
  107.     x = a;
  108.     x.a++;
  109.     return x;
  110. }
  111.  
  112.  
  113. define surd_dec(a)
  114. {
  115.     local    x;
  116.  
  117.     x = a;
  118.     x.a--;
  119.     return x;
  120. }
  121.  
  122.  
  123. define surd_neg(a)
  124. {
  125.     local obj surd    x;
  126.  
  127.     x.a = -a.a;
  128.     x.b = -a.b;
  129.     return x;
  130. }
  131.  
  132.  
  133. define surd_mul(a, b)
  134. {
  135.     local obj surd    x;
  136.  
  137.     if (!istype(b, x)) {
  138.         x.a = a.a * b;
  139.         x.b = a.b * b;
  140.     } else if (!istype(a, x)) {
  141.         x.a = b.a * a;
  142.         x.b = b.b * a;
  143.     } else {
  144.         x.a = a.a * b.a + surd_type * a.b * b.b;
  145.         x.b = a.a * b.b + a.b * b.a;
  146.     }
  147.     if (x.b)
  148.         return x;
  149.     return x.a;
  150. }
  151.  
  152.  
  153. define surd_square(a)
  154. {
  155.     local obj surd    x;
  156.  
  157.     x.a = a.a^2 + a.b^2 * surd_type;
  158.     x.b = a.a * a.b * 2;
  159.     if (x.b)
  160.         return x;
  161.     return x.a;
  162. }
  163.  
  164.  
  165. define surd_scale(a, b)
  166. {
  167.     local obj surd    x;
  168.  
  169.     x.a = scale(a.a, b);
  170.     x.b = scale(a.b, b);
  171.     return x;
  172. }
  173.  
  174.  
  175. define surd_shift(a, b)
  176. {
  177.     local obj surd    x;
  178.  
  179.     x.a = a.a << b;
  180.     x.b = a.b << b;
  181.     if (x.b)
  182.         return x;
  183.     return x.a;
  184. }
  185.  
  186.  
  187. define surd_div(a, b)
  188. {
  189.     local x, y;
  190.  
  191.     if ((a == 0) && b)
  192.         return 0;
  193.     obj surd x;
  194.     if (!istype(b, x)) {
  195.         x.a = a.a / b;
  196.         x.b = a.b / b;
  197.         return x;
  198.     }
  199.     y = b;
  200.     y.b = -b.b;
  201.     return (a * y) / (b.a^2 - surd_type * b.b^2);
  202. }
  203.  
  204.  
  205. define surd_inv(a)
  206. {
  207.     return 1 / a;
  208. }
  209.  
  210.  
  211. define surd_sgn(a)
  212. {
  213.     if (surd_type < 0)
  214.         quit "Taking sign of complex surd";
  215.     if (a.a == 0)
  216.         return sgn(a.b);
  217.     if (a.b == 0)
  218.         return sgn(a.a);
  219.     if ((a.a > 0) && (a.b > 0))
  220.         return 1;
  221.     if ((a.a < 0) && (a.b < 0))
  222.         return -1;
  223.     return sgn(a.a^2 - a.b^2 * surd_type) * sgn(a.a);
  224. }
  225.  
  226.  
  227. define surd_cmp(a, b)
  228. {
  229.     if (!istype(a, surd__))
  230.         return ((b.b != 0) || (a != b.a));
  231.     if (!istype(b, surd__))
  232.         return ((a.b != 0) || (b != a.a));
  233.     return ((a.a != b.a) || (a.b != b.b));
  234. }
  235.  
  236.  
  237. define surd_rel(a, b)
  238. {
  239.     local x, y;
  240.  
  241.     if (surd_type < 0)
  242.         quit "Relative comparison of complex surds";
  243.     if (!istype(a, surd__)) {
  244.         x = a - b.a;
  245.         y = -b.b;
  246.     } else if (!istype(b, surd__)) {
  247.         x = a.a - b;
  248.         y = a.b;
  249.     } else {
  250.         x = a.a - b.a;
  251.         y = a.b - b.b;
  252.     }
  253.     if (y == 0)
  254.         return sgn(x);
  255.     if (x == 0)
  256.         return sgn(y);
  257.     if ((x < 0) && (y < 0))
  258.         return -1;
  259.     if ((x > 0) && (y > 0))
  260.         return 1;
  261.     return sgn(x^2 - y^2 * surd_type) * sgn(x);
  262. }
  263.  
  264. global lib_debug;
  265. if (lib_debug >= 0) {
  266.     print "obj surd {a, b} defined";
  267.     print "surd(a, b) defined";
  268.     print "surd_print(a) defined";
  269.     print "surd_conj(a) defined";
  270.     print "surd_norm(a) defined";
  271.     print "surd_value(a, xepsilon) defined";
  272.     print "surd_add(a, b) defined";
  273.     print "surd_sub(a, b) defined";
  274.     print "surd_inc(a) defined";
  275.     print "surd_dec(a) defined";
  276.     print "surd_neg(a) defined";
  277.     print "surd_mul(a, b) defined";
  278.     print "surd_square(a) defined";
  279.     print "surd_scale(a, b) defined";
  280.     print "surd_shift(a, b) defined";
  281.     print "surd_div(a, b) defined";
  282.     print "surd_inv(a) defined";
  283.     print "surd_sgn(a) defined";
  284.     print "surd_cmp(a, b) defined";
  285.     print "surd_rel(a, b) defined";
  286.     print "surd_type defined";
  287.     print "set surd_type as needed";
  288. }
  289.